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Abstract 

Using Monte Carlo simulations we study two-dimensional prey-predator systems. Measuring the 
variance of densities of prey and predators on the triangular lattice and on the lattice with eight 
neighbours, we conclude that temporal oscillations of these densities vanish in the thermodynamic 
limit. This result suggests that such oscillations do not exist in two-dimensional models, at least 
when driven by local dynamics. Depending on the control parameter, the model could be either in 
an active or in an absorbing phase, which are separated by the critical point. The critical behaviour 
of this model is studied using the dynamical Monte Carlo method. This model has two dynamically 
nonsymmetric absorbing states. In principle both absorbing states can be used for the analysis of 
the critical point. However, dynamical simulations which start from the unstable absorbing state 
suffer from metastable-like effects, which sometimes renders the method inefficient. 



I. INTRODUCTION 



Nonequilibrium, many-body systems were usually regarded as one of the classical domains 
of physics. However, recently is is becoming clear that such systems are of interest also in 
other disciplines such as biology economy 0, or sociology ||. Consequently, notions of 
cooperative effects, formation of order or phase transition, which were typically restricted 
to physics, proliferate other scientific disciplines. 

An example of multidisciplinary problem, which has application both in biological pop- 
ulation dynamics and in hetero- or homogeneous catalysis, is the existence of temporarily 
periodic solutions in the so-called prey-predator systems. A closely related problem is syn- 
chronization, which appears in various biological contexts such as heart beats, wake-sleep 
cycles and many other Q. Earliest mathematical model of the temporal oscillations in inter- 
acting populations was proposed by Lotka and Volterra, who described populations of preys 
and predators in terms of nonlinear differential equations || . However, description of inter- 
acting populations in terms of ordinary differential equations, by necessity, introduces con- 
siderable simplifications. In particular it neglects fluctuations and spatial inhomogeneities 
which are certainly present in such systems. One possibility to take such effects into account 
is to describe such systems using lattice models. The first approach of this kind was made 
by Chate and Manneville ||. Subsequently, more complicated versions of their model were 
studied 0. Although in some cases these studies confirmed existence of periodic in time 
solutions, these models were driven by synchronous dynamics. More realistic models with 
asynchronous dynamics were also examined [|S], |IU|| . As we already mentioned, closely 
related models are studied in the context of catalysis (Tl| . 



On general grounds one expects that fluctuations in low- dimensional systems are strong 
enough to destroy such temporal oscillations in the thermodynamic limit. An important 
question is: what is the critical dimension d c above which such oscillations will survive in 



the thermodynamic limit. Certain arguments were given by Grinstein et al. that d c = 2 [12 



Numerical results for both synchronous || [13| and asynchronous models [14 seem to confirm 
that temporal oscillations exist but only for d > 2. Recently, however, numerical simulations 
by Rosenfeld et al. suggest that in a model with smart preys and predators oscillatory 
behaviour might exist even in d — 2 systems [[OJ. They argued that the oscillatory phase 
is induced by a certain dynamical percolation transition. If this percolative mechanism 
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would be of more general validity, then the oscillatory phase should exist also in other two- 
dimensional models. One would then argue that the fact that this phase was not observed 
in some earlier simulations |14| was due to for example too small range of interactions of 
preys and predators. 

In the present paper we examine two-dimensional versions of a prey-predator model |TJj 
for the increased range of interaction. However, our simulations show that even with the 
increased range of interaction, an amplitude of oscillations vanishes in the thermodynamic 
limit. Let us notice that in three-dimensional version of the model, oscillations exists even 



for the nearest-neighbour interactions Hl4fl . Thus, our result strongly indicates that, at least 
for our model, there are no oscillations in the case d = 2. 

Usually lattice prey-predator models have an absorbing state. Consequently, they might 
undergo a phase transition in the steady state between active and absorbing phases, which is 



expected to belong to the directed percolation universality class fllEf . Our model is charac- 
terized by two absorbing states. However, these states are asymmetric and only one of them 
is typically reached by the model's dynamics. The second absorbing state can be considered 
as dynamically unstable. Because of this asymmetry the model behaves as if it has only one 
absorbing state and thus should exhibit directed percolation criticality, which was confirmed 
for the one-dimensional case | I7| . Critical behaviour of models with absorbing states can be 
studied using the so-called dynamical Monte Carlo method |H5|]. In this method we prepare 
the system in an absorbing state and locally activate it. Statistics of such runs provides 
a very efficient method to study the critical behaviour of such models. This method was 
already applied to a number of models with a single absorbing state or with double but 
symmetric ones ||20|| . We are, however, not aware of studies for models with asymmetric 
absorbing states. We show that when dynamical simulations start from an unstable absorb- 
ing state they suffer from metastable-like effects, which in the two-dimensional case renders 
the method inefficient. In the one- dimensional case this effect is much weaker and both 
absorbing states can be used to extract information about the critical point. 

In section II we define the model, briefly mention the results of the simple mean-filed ap- 
proximations, and study the fluctuations of densities of populations. In section III we present 
the results of the dynamical Monte Carlo method. Section IV contains our conclusions. 
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II. MODEL AND ITS STEADY-STATE SIMULATIONS 



In our model on each site of a ci-dimensional cartesian lattice of linear size N we have 
a four-state variable = 0,1,2,3, which corresponds to the site being empty (e = 0), 
occupied by a prey (e = 1), occupied by a predator (e = 2) or occupied by a prey and a 



predator (e = 3). Dynamics of this model is specified as follows [14]: 

(i) Choose a site at random (say z-th). 

(ii) With the probability r (0 < r < 1) update a prey at the chosen site, provided that 
there is one (i.e., e = 1 or 3); otherwise do nothing. Provided that at least one neighbour 
of the chosen site is not occupied by a prey (ie., e = or 2), the prey (which is to be 
updated) produces one offspring and places it on the empty neighbouring site (if there are 
more empty sites, one of them is chosen randomly). Otherwise (i.e., when there is a prey 
on each neighbouring site) the prey does not breed (due to overcrowding). 

(iii) With the probability 1 — r update a predator at the chosen site, provided that there 
is one (i.e., e = 2 or 3). Provided that the chosen site is occupied by a predator but is not 
occupied by a prey (e = 2), the predator dies (of hunger). If there is a prey on that site 
(i.e., e = 3), the predator survives and consumes the prey from the site it occupies. If there 
is at least one neighbouring site which is not occupied by a predator, the predator produces 
one offspring and places it on the empty site (chosen randomly when there are more such 
sites). 



To complete the description of this model we have to specify what are the neighbouring 
sites, i.e., sites where offsprings can be placed. In previous studies of this model these were 
just nearest neighbours fl4] , [17| . In the present paper of our main concern are d = 2 models 
with further (but finite) range. 

Steady-state description of our model is given in terms of densities of preys x and preda- 
tors y defined as 

x = J7d XX<W + S a,a), ^^Efe + M- (!) 

i i 

where summation is over all N d sites i and 5 is Kronecker's 5-function. From the above 
rules it follows that the model has two absorbing states. The first one (PREY) is filled with 
preys (x — 1, y — 0) and the second one (EMPTY) is empty (x — y — 0). For large enough 
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r, both populations coexist and the model is in the active phase (x > 0, y > 0). When 
the update rate of preys r decreases, their number become to small to support predators. 
For sufficiently small r predators die out and the system quickly reaches the absorbing state 
PREY. In the simplest case, when the neighbouring sites are nearest neighbours, the phase 
transition between active and absorbing phase was observed at positive r for d = 1, 2, 3 |T7] . 
For d = 1 (linear chain) Monte Carlo simulations show that, as expected, the phase transition 
belongs to the directed percolation universality class. 

To get additional insight into the behaviour of the model we can use a mean-field ap- 
proximation. In its simplest version one describes the model solely in terms of densities 
x and y (which are strictly speaking time averages quantities defined in Eq. (fj)). From 
the dynamics of the model and typical mean-field assumptions one can easily arrive at the 



following equations [14, 17 



dx 

— = rx{l-x w )-(l-r)xy, (2) 

J = (1 - r)xy{\ - y w ) - (1 - r)y(l - x). (3) 

where w is the number of neighbouring sites (see the dynamical rule (ii) and (iii)). Pre- 
dictions of such an approximation is, however, substantially different from the behaviour of 
the model as observed in Monte Carlo simulations [Q. In particular, approximation (0)-(0) 
fails to predict a phase transition between active and absorbing phases at positive r and for 
any dimension. Moreover, there is no indication of the oscillatory phase as observed in the 
d = 3 case. 

An improved version of the mean-field approximation can be obtained introducing a third 
variable z, which denotes a density of sites occupied by a prey and a predator (e = 3). Then, 
the mean-field equations are written as: 

— = rx{l-x w )-(l-r)z, (4) 

(lb 

f t =(l-r)z(l-y w )-(l-r)(y-z), (5) 

dz _ rxjl - x w )(y - z) (1 - r)z{l + z - x - y){l - y w ) , w 
at 1 — x 1 — y 

In the approximation (@)-(@) the density of sites that are occupied by a prey and a predator is 
simply given by the product xy. In the approximation (§)-(§[) this is an independent variable, 
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FIG. 1: The variance of the density of preys as a function of 1/N for the model on the triangular 
lattice. Simulations were made for r = 0.4(D), 0.3(*),0.2(x), and 0.15 (+). 

whose time evolution follows from the dynamical rules of the model. Numerical solution of 



(ED~(§) shows that this approximation is indeed more accurate In particular, it predicts 
that for w > 4 (which could be interpreted as d > 2) there is a range of r in the active phase 
where nonvanishing oscillations exist. For w = 4 (i.e., square lattice) this is still at odd with 
Monte Carlo simulations [0. Let us emphasize that both Monte Carlo simulations and 
mean-field calculations are essentially independent on the initial configuration (provided, of 
course, we do not start from an absorbing state). 

To examine the oscillatory behaviour, one can measure the standard deviation of the 
densities defined as a T - 



< (x— < x >) 2 > and analogously for the densities of predators 
y. The symbol < ... >, denotes time average in the steady state. It was found that for 
N — > oo a = a x always decreases to zero for d — 1,2 fH]. However, for d = 3 there is a 



certain range of r in the active phase and such that a remains finite for TV — > oo. Such a 
behaviour indicates that finite- amplitude oscillations survive in the thermodynamic limit. 

To check, whether an increased range of interactions might stabilize oscillations in two- 
dimensional systems, we simulated our model on the triangular lattice and on the lattice with 
8 neighbours (square lattice with nearest- and next-nearest neighbours). On such lattices 
the parent site that is supposed to breed looks for an empty site among 6 or 8 neighbouring 
sites, respectively. On such lattices we also observed a phase transition into absorbing phase 
around r ~ 0.06 — 0.08. Similarly to other lattices fl4[], for r larger but close to the transition 
point fluctuations of densities exhibit oscillatory-like behaviour. However, analysis of the 
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FIG. 2: The variance of the density of preys as a function of 1/N for the model on the lattice with 
nearest- and next-nearest neighbours. Simulations were made for r = 0.35(*), 0.3(x), and 0.25 
(+)• 

variance of this fluctuations shows (see Fig. p] and Fig. |2|) that in the thermodynamic limit 
amplitude of oscillations vanishes. Simulations were made for N < 1000 and simulation 
time was long enough to ensure that error bars are smaller than the plotted symbols. 

Let us notice that for the lattice with 8 neighbours the coordination number is greater 
that in the three-dimensional (cubic lattice) case. While in the former case oscillations die 



out in the thermodynamic limit, they survive in the latter one [14]]. It indicates that in 
two-dimensional systems oscillations are unlikely. However, such a conclusion, if true, is 
applicable only to models equipped with local dynamics. 

One factor which is omitted in our approach is diffusion. Our organisms diffuse only 
effectively through the breeding process. In related two-dimensional models, which were 



studied in the context of heterogeneous catalysis, diffusion was taken into account pi] . It 
seems, however, that as long as the diffusion constant is finite, the amplitude of oscillations 
vanishes in the thermodynamic limit. 



III. DYNAMICAL MONTE CARLO 



The essence of this method is to prepare the system in an absorbing state and to activate 
it locally [[19]. Quantities of the main interest are the probability P(t) that the activity did 
not die out until time t (the unit of time is defined as a single, on average, update of each 
site) and the average number of active sites N(t) (averaged over all runs). One expects, that 
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FIG. 3: The average number of predators N(t) as a function of time t for the model on the square 
lattice with the absorbing state PREY. Simulations were made for (from top) r = 0.1099, 0.1098, 
0.1097 and 0.1096. The dotted line has a slope 0.23 which corresponds to the exponent r\ in (2+1) 
DP HJ. 
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FIG. 4: The survival probability of predators P(t) as a function of time t for the model on the 
square lattice with the absorbing state PREY. Simulations were made for (from top) r = 0.1099, 
0.1098, 0.1097 and 0.1096. The dotted line has a slope 0.451 which corresponds to the exponent 5 
in (2+1) DP @. 

at criticality these quantities have the power-law behaviour: P(t) ~ t~ s , N(T) ~ t v , where 
5, and rj are critical exponents characteristic to a given universality class. Off the critical 
point P(t) and N(t) deviate from the power-law behaviour. 

We applied the dynamical Monte Carlo method to the two-dimensional model with near- 
est neighbour interactions (square lattice). In this case the steady-state calculations show 
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that the model undergoes a transition at r = r c ~ 0.11. Our results in the case of the ab- 
sorbing state PREY are shown in Fig. 0-Fig. |j. As an initial configuration we put all sites in 
a prey state (e = 1), except a one, randomly selected, site which is in a prey-predator state 
(e = 3). Actually, in our model there are two quantities which scale approximately in the 
same way, namely, the number of predators and the number of active sites (a site is active if 
it contains a predator, a prey and predator or a prey which is surrounded by at least one site 
without a prey). To calculate P(t) and N(t) we considered only the number of predators 
(when this number becomes zero we stop a given trial) but qualitatively the same results 
are obtained when the number of active sites is considered. Using dynamical Monte Carlo 
method it is important to ensure that the spreading activity do not reach the boundary of 
the system. We checked that the value which we used (N = 2000) was sufficiently large. 
From the behaviour of P(t) and N(t) we conclude that the phase transition is located at 
r = r c = 0.10975(5), which is far more accurate than the steady-state estimation. Moreover, 
one can see that our data upon approaching the critical point have a power-law behaviour 
with exponents consistent with (2+1) directed percolation [ f2~0[ . 



Now, let us apply the dynamical Monte Carlo in the case of the EMPTY absorbing state. 
As an initial configuration we put all sites in an empty state (e = 0), except a one, randomly 
selected, site which is in a prey-predator state (e = 3). Numerical results are shown in 
Fig. |5|-Fig. H In these simulations the system is in the active phase (r > r c ) and relatively 
far from the critical point. Although the survival probability P(t) saturates for large t its 
asymptotic values are very low (Fig. ||). For example, to confirm with this method that 
for r = 0.12 the system is in the active phase tens of millions of runs must be made. The 
reason for that is that for most of the runs predators die out very quickly. Only very few 
runs generate sufficiently many preys which can support the population of predators. If 
the population of predators survives for a certain time, then, most likely, it will survive for 
an infinitely long time spreading through the system (see Fig. ||). We expect that such a 
behaviour exists for any r > r c . However, for r close to r c the long-time survival events are 
extremely unlikely, which renders this method inefficient. 

Our model in this case resembles a metastable system |]2"2"fl . Such a system can exit the 
metastable state only through formation of a sufficiently large seed of the stable phase. In 
our case the EMPTY absorbing state corresponds for r > r c to a metastable state. Only 
when a sufficiently large island of the active state occurs, the system is irreversibly driven 
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FIG. 5: The survival probability of predators P(t) as a function of time t for the model on the 
square lattice with the absorbing state EMPTY. Simulations were made for (from top) r = 0.25, 
0.2, 0.15 and 0.12. 

toward the active (stable) state. 

However, this metastable effect is much weaker in one-dimensional systems. In this case 
the critical point is at r = r c = 0.491(2) |I7J. Dynamical Monte Carlo simulations for both 
PREY and EMPTY absorbing states are relatively efficient (see Fig. In the case of the 
EMPTY absorbing state only a mild decrease of N(t) is seen for small t. 

We do not fully understand why metastability is much weaker in the d = 1 case. Most 
likely it is due to a much easier formation of a sufficiently large seed of a stable phase. 
Outward spreading of such an active seed resembles a biased random walk and thus grows 
steadily in time. In the d = 2 case, the formation and growth of the active seed is more 
complicated. 

IV. CONCLUSIONS 

In the present paper we studied oscillatory and dynamical properties of the two- 
dimensional prey-predator systems. Our results show that even for the increased range 
of interactions an amplitude of oscillations vanishes in the thermodynamic limit. It might 
suggests, as found by Grinstein et al. in a related context [0 , that oscillations might persist 
but only for d > 2. 

In addition to that, we examined the validity and efficiency of the dynamical Monte Carlo 
method for a model with two asymmetric absorbing states. It turns out that for an unstable 
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FIG. 6: The average number of predators N(t) as a function of time t for the model on the square 
lattice with the absorbing state EMPTY. Simulations were made for (from top) r = 0.25, 0.2, 0.15 
and 0.12. 
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FIG. 7: The average number of predators N(t) as a function of time t for the model on the 
one-dimensional lattice with the absorbing state PREY and EMPTY. For each absorbing state 
simulations were made for (from top) r = 0.493, 0.491 and 0.489. The slope of the critical (central) 
lines is very close to the (1+1) directed percolation value r\ = 0.314 f2Cfl . 

absorbing state in the two-dimensional system this method is very inefficient. On the other 
hand, both absorbing states can be used within this method in a one-dimensional model. 
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